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ABSTRACT 

We propose a monitoring indicator of the normality of the out- 
put of a gravitational wave detector. This indicator is based 
on the estimation of the kurtosis (i.e., the 4th order statisti- 
cal moment normalized by the variance squared) of the data 
selected in a time sliding window. We show how a low cost 
(because recursive) implementation of such estimation is pos- 
sible and we illustrate the validity of the presented approach 
with a few examples using simulated random noises. 

Four large-scale detectors ^ of gravitational waves (GWs) 
are on the point to take their first scientific data. They all rely 
on the principle used in the Michelson experiment: measure 
the relative length 8L/L of the two perpendicular arms (each 
formed by two suspended masses) of the detector. The goal 
is to reach a measurement sensitivity (i.e., decrease the noise 
level) such that the small changes of 5L caused by the GWs 
emitted from astrophysical sources (such as the coalescing bi- 
naries of neutron stars or black holes) can be detected when 
they pass through the instrument. 

Since only a small number of such events are likely to be 
observed, the problem consists from the data analysis view- 
point, in looking for rare transients appearing in the detector 
output. For the coalescing binaries mentioned above, this will 
be done by implementing a bank of matched filters: each fil- 
ter correlate the data with the expected GW emitted from a 
binary and we repeat this operation for a large number of pos- 
sible targets. Assuming the template waveform are reliable, 
the matched filtering approach can be shown to be optimal in 
the Neymann-Pearson sense provided that the noise is Gaus- 
sian. Therefore, it is crucial to monitor the noise Gaussianity 
and locate any departures from this nominal hypothesis. 

There are many ways for the noise to departs from Gaus- 
sianity, however not all of them are relevant for the problem 
of GW detection. One of the possibilities is to have a noise 
probability density function (PDF) with heavier tails than the 
Gaussian bell curve. This particular discrepancy is a problem 
because it causes a increase of the false detection rate (i.e., the 
large values in the tails make the matched filter triggers more 
often). 

The kurtosis [||| is a well-known measurement of the de- 
cay rate of the PDF in the tails. This motivates us to use 



the kurtosis as an index measuring normality. We define the 
mean of the signal x(t) by jtzi(t) = E[x(t)] where IE[-] is 
the expectation operator. The central moments [^] of x(t) are 
given by fi r (t) = E[(x(t) — n 1 (t)) r ] where r > 2 is the or- 
der. The kurtosis is defined by the following ratio K^t) = 

Because the analysis must be done in real-time, the nor- 
mality test should not be computationally expensive. We pro- 
pose here an efficient (because recursive) implementation of 
the kurtosis estimation in a short-term and sliding observa- 
tion window. Another recursive estimator of the kurtosis was 
proposed in [||] and required to have zero-mean signals. The 
presented approach here works also with non-centered sig- 
nals. 

The outline of the paper is as follows. We choose a sim- 
ple mathematical structure for the short-term and recursive 
estimation of the central statistical moments. In Sect, [lj we 
identify in the selected family which estimators of /ii, [i 2 and 
/i4 have a vanishing bias. We then show in Sect. || that an 
adequate Taylor approximation of the ratio of the unbiased 
estimates of \i± and fi 2 obtained previously yields the recur- 
sive estimate of K4 and we detail its computation algorithm. 
Finally, we apply in Sect. || the proposed estimator to a few 
illustrative cases and we explain how it is used in practice for 
the monitoring of the normality of the GW detector output. 

1. RECURSIVE ESTIMATES WITH VANISHING 
BIAS 

1.1. Recursive estimate of the mean and variance 

We assume that the signal x(t), t 6 Z (using a unit sampling 
rate) is locally stationary (i.e., its statistical moments do not 
change during a finite time period 5t sta t) and ergodic (i.e., its 
statistical moments can be estimated from its samples). The 
mean of x(t) can be estimated with the following weighted 
average of the data selected by the window wi (t) 



fa(t) =Ciy]wi(t- n) x(n), 

n=0 



(1) 
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where the duration of wi(t) is smaller than 8t s tat- If the win- 
dow in use is of exponential type, i.e. w\(t) = a\, this esti- 
mator can be equivalently calculated recursively with 

(l 1 {t)=aifl 1 (t-l) + C 1 x(t). (2) 



The problem is to find the constants a\ > and C\ such 
that (ii(t) is asymptotically unbiased i.e., E[/2i(i)] = /ii 
when t — > +00. The bias can be calculated directly from 
the definition of the estimator, yielding 



E[AiW] = Cimi 



1 



,t+l 



1 — a\ 



1 — ai 



/ii when t — ► +00 



and assuming ai < 1. We conclude that (ii(t) is an unbiased 
estimator of the mean if C\ = 1 — a%. 

The same method can be applied to find recursive and un- 
biased estimators of the higher order moments, based on the 
following choice of expression : 

t 

fi r {t) = C r w r (t — n) (x(n) — fii(n — l)) r , (3) 

71=0 

where r > 2. Similarly to the mean, a recursive implementa- 
tion is possible if we choose w r (t) = a*. We restrict to the 
interval < a r < 1. 

When r = 2, we essentially average the squared differ- 
ences to the estimated mean in a sliding time window defined 
by W2(t). We evaluate the bias of (12(f) by taking the expec- 
tation of (|3|). We first evaluate : 



E[(>(n)-Ai(n-1)) 




and setting Ci = 1 — a\ (i.e., set the bias of fix to zero), the 
summation leads to : 



EM*)] = 



2C 2 ^2 I — a. 



t+i 
2 



1 + ai 1 — a 2 



C2 Mi - 1 , _ M2 



„t+i 



2(t+l) 



1 + ai 
2C 2 



a 2 - oj 



(l + ai)(l-aa) 



M2, 



when £ — > +00 provided that a-y and a 2 < 1. Consequently, 
the bias of /t 2 (t) is zero when C2 = (1 + &i)(l — a 2 )/2. 

1.2. Extension to the 4th order 

We proceed to the fourth order as previously, starting from the 
definition (||) with r = 4. The evaluation of the bias of jii(t) 
requires the calculation of 

E[(x(n)-fi 1 (n-l)) i ] = (l + f 4 ) f i 4 
+ 4(1 -h-h + /i/ 3 )MiM3 + 3(2/ 2 + fl - 
+ 6(1 - 2/ x + h + fl - 2A/ 2 + A 2 / 2 )m?M2 

+ (1 - 4.A + 6A 2 - Afl + U)& 

where we defined f k = C\ (1 — aj n )/(l — af). 

Setting C\ = 1 — ai so that the estimate of the mean is 
not biased, and making the summation for n = . . . t, we get 



E[£ 4 (*)] 



C 4 



1 — a 4 



(fc 4 ^4 + k 2 ^t) + r(t). 



(4) 



where r{t) is a complicated sum of integer powers of terms of 

the form Ka^ and Ka^ +1 ^ where K € K.. The constants 
can be expressed as 



2(1 -01 + 2a 2 ) 
(l + ai)(l + a 2 ) 



6(l-ai)(l + a 1 + 2a 2 ) 
(l + ai ) 2 (l + a 2 ) ■ 



It is interesting to note that the first term in (Q) does not 
depend on ^3. 

If 01 and a 4 < 1, the function r(t) goes to when t tends 
to +00 so that the expectation of fi^t) depends in a simple 
manner of /14 (the objective value) and Assuming that /i 2 
is known, we can set the bias of ji&(t) to a simple constant 
offset if we choose the coefficient of \i± in (Q) equal to 1, i.e. 
we set C4 = (1 — a 4 )/fc4 in the following. 

2. TAYLOR APPROXIMATION OF THE KURTOSIS 
ESTIMATOR 

The results of the previous Section motivate us to propose 
the following estimator ki(t) of the kurtosis, obtained by di- 
viding (Q) by /1 2 and replacing the variance by its recursive 
estimate : 

k 2 

K 4 (t) = K±(t) - — , (5) 
K4 

where Ri(t) = /i 4 (i) / p%(t) is a heuristic estimator which we 
correct in (|5|) by suppressing the offset. 

We choose the same window for all the involved estima- 
tors (i.e., we set 04 = a 2 = ax), so that 



K4,(t) 



ailHjt - 1) + Cijxjt) - /Ti(t - l)) 4 
(a lf i 2 (t - 1) + dixit) - Mt - wy 



For sufficiently large duration of the window (ai — * 1), 
we can treat C\ as an epsilon and approximate this ratio to the 
first two terms of its Taylor expansion (this idea was inspired 
by [f|]X which leads to the following expression : 

« 4 (f) = (1 + Ci - 2Ci5 2 x{t)) Ki{t - 1) 

+ C 1 8 A x{t) + 0(Cl), (6) 

where 5 2 x(t) = (x(t) — /i 1 (t — l)) 2 //i 2 (t — 1) computes a 
normalized distance of the current signal sample to the mean. 
The correct estimate of the kurtosis proposed in (|5|) is ob- 
tained by subtracting the offset given by 



h _ -3Ci(4- 5Ci +2C 2 ) 
h ~ (2-Ci)(2-3Ci +2C* 2 ) 



= -3d + 0(C 2 ), 



to the approximation in (g). Note that since C\ -C 1, the 
offset can be neglected in most practical situations because 
«4 > 1. 

The recursive estimation of the kurtosis obtained in (^) is 
intuitively appealing : if we replace the estimators of the mean 
and variance by their actual values in the definition of 5 2 x(t), 
then noting that E[d 2 x(t)] = 1 and IE[c5 4 x(£)] = Ki(t), we 
can write 



K4(t) Ri Cli Ki(t — 1) + CiKi(t). 



(7) 



Using the equivalence between ([[]) and (||), we conclude 
that eq. (Q) may be seen as a local average of the crude esti- 
mation of «4(i) made by 8 A x(t). 

Eq. (|6j) yields a recursive estimation scheme which is de- 
scribed by the pseudo-code in Tab. [j] 

3. VALIDATION AND PRACTICAL USE 

In this section, we make various numerical checks of the pro- 
posed method. In all the simulations, the estimator ki(t) as 
defined in eqs. (^) and (^|) is computed with C\ = 2.9912 x 
10~ 3 . This corresponds to a window duration of about 20 s if 
a sampling rate f s = 50 Hz is assumed. 
Check #1: what is the bias of the estimator? — We answer 
this question in the case of a Gaussian noise for which n± = 3. 
Figure [l] (top) shows the histogram of «4(t) computed with 
a simulated (zero-mean, unit variance and white) Gaussian 
noise. With this histogram, we evaluate the expectation of 
K4(t) to be equal to 3.0006 with a bin size of ±0.03. 
Check #2: is the estimator useful for detecting noises with 
heavy tails? — Figure ^ presents the result of the estima- 
tion of the kurtosis of an evolving mixture of Gaussian and 
Laplacian noises. The kurtosis of a Laplace random variable 
is equal to K4 = 6 indicating that the tails of this distribu- 
tion oc exp(— \/2|x|) [@] decrease slowly as compared to the 
Gaussian ones. The two noises are linearly combined ; the 
weight coefficient of the Laplacian noise increases from to 
1 (and reverse for the Gaussian noise) following a linear func- 
tion of time. An excess of kurtosis (i.e., £4(7;) > 4) appears 
starting from 1 1=3 80s (time at which the Laplacian noise starts 
to dominate) showing that we can answer positively to the 
question. 

Check #3: effect of non-stationarities — Figure |]illustrates 
how the recursive estimator of the kurtosis performs with a 
simulated Gaussian noise of changing mean and variance. Af- 
ter a transient period roughly equal to the window length, 
£4(7;) tends to the correct value which is 3. Each change of 
the mean or variance is seen as a non-Gaussianity (i.e., large 
values of the kurtosis). The reason is that the hypothesis of 
local stationarity required for a correct estimation (see Sect. 
|]) is not satisfied at the discontinuity points. 
Practical use with gravitational wave data — For technical 
reasons related to the common data format used by the gravi- 
tational wave detectors, it is convenient to fix the output rate 
of the monitoring indicators to one sample per (GPS) second 
of data (also referred to as frame). This applies to the normal- 
ity index we would like to set up. 

Since the variance of k^(t) is difficult to obtained, we can- 
not compute a confidence interval which would be required 
to conclude on the normality of the data from the estimator 
value. We remedy this with the following scheme : 

1. choose arbitrarily a threshold 7/ (e.g., r\ = 4), 

2. associate a warning flag to each frame of data where an 
excess of kurtosis (i.e., £4 > rj) has been observed at 
least once, 



3. compute the rate of triggered frames (over periods cor- 
responding to the typical duration of a GW as seen by 
the detector), 

4. compare this rate to the one evaluated with Monte-Carlo 
simulations using a Gaussian noise with similar spec- 
tral characteristics than the signal being observed. 

Figure [l] (bottom) gives the expected value (with error 
bars) of this rate if the signal is Gaussian and white, for thresh- 
olds taken between 2.4 and 4. For instance, if we fix 77 = 4, 
then rates of triggered frames larger than 1 % indicate the pres- 
ence of a heavy-tailed noise in the data. 



Table 1. Pseudo-code for the computation of k 4 (t). The 

algorithm requires a total number of 16 floating point opera- 
tions (10 multiplications and 6 additions) to compute the next 
kurtosis value from the previous one. The memory usage is 
restricted to three registers of real numbers. As a comparison, 
a sliding estimate using k-Statistics ^ would need at least 
a register of the same size than the window. Note that there 
are several possible initializations of the three registers in line 
5. In our simulations, we set them to 0, 1 and respectively. 
This initialization affects essentially the transient period at the 
beginning of a computation. 



choose a window duration -> define al 

Cl=l-al 

C2= (l-al~2) 12 

bias=-3 C_l 

init mul_last, mu2_last, k4_bar_last 

while (data available) 

x = get next data sample 

mul= al mul_last + CI x 
dx2= (x-mul_last ) "2 
mu2= al mu2_last + C2 dx2 
dx2= dx2/mu2_last 

k4_bar= (1+C1-2 CI dx2) k4_bar_last ... 

... + CI dx2"2 
kappa4= k4_bar + bias 

send kappa4 to output 

mul_last= mul 
mu2_last= mu2 
k4_bar_last=k4_bar 
end 
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Fig. 1. Probability density function of k± and rate of trig- 
gered frames computed with Gaussian noise |@]. The re- 
sults presented in this figure were obtained using 50 streams 
of simulated (zero mean, unit variance and white) Gaussian 
noise. Each data stream contains 30,000 samples (i.e., 600s 
if the sampling rate is f s = 50 Hz), top: the empiri- 
cal histogram presented in this diagram gives an estimation 
(bounded by error bars) of the PDF of k 4 (t) (see Sect. [5] 
for details), bottom: for a threshold -q taking values between 
2.4 and 4, we present the percentage of frames of data where 
ki(t) > j] at least once. To make this computation, the data 
was divided into 600 chunks of 1 second duration (each of 
them defining a frame) and the first 20 frames (i.e., equivalent 
to the window duration) of each trial were removed. 
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Fig. 2. Recursive estimate of the kurtosis of an evolving 
mixture of Gaussian and Laplacian noises top: the test 
signal we use is a mixture of noises (both of unit mean and 
variance) having resp. a normal PDF A/"(l, 1) and a Laplace 
PDF oc exp(— \/2\x — 1|). The two noises are combined lin- 
early. The weight coefficient of the Laplacian noise increases 
from to 1 (and reverse for Gaussian) following a linear func- 
tion of time, bottom: we see in this plot that an excess of 
kurtosis ki{t) (see Sect. [3] for details) appears starting from 
t w 80s (time at which the Laplacian noise starts to domi- 
nate). 
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Fig. 3. Recursive estimate of the kurtosis of a Gaussian 
noise of changing mean and variance [|[. top: this is the 
same simulation as in Fig. || with a white Gaussian noise 
whose mean and variance are set resp. to the following values 
(1, -2, 1) and (1, 1,4) during the three periods (t < 54 s, 
54 < t < 108 s, t > 108 s). bottom: after a transient period 
roughly equal to the window length, k^lt) tends to the correct 
value which is 3. Each change of the mean or variance is seen 
as a non-Gaussianity (i.e., large values of the kurtosis). 



